package fem2.element;

import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Node;
import fem2.element.crack_sda.SDAElement;
import fem2.enu.MeshPartType;
import fem2.material.cohesive_law.CohesiveLaw;
import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.math.LinAlg;

import java.util.ArrayList;

import math2.MathUtilities;

/**
 * constant jump element based on Linder-Armero paper
 * 
 * @author hbui
 * 
 */
public class ConstantJumpLAElement extends SDAElement {

	private IVector Xi;

	private double crackLength;
	private double elementArea;
	private double[] normal;
	private double[] tangent;

	private ArrayList<Integer> posDomain;
	private ArrayList<Integer> negDomain;

	public ConstantJumpLAElement(MeshPart mp, MaterialModel mm, CohesiveLaw cohesiveLaw) {
		super(mp, mm, cohesiveLaw);
		setCrackMeshPartType(MeshPartType.LINE2);
		setIntegrationOrderOnCrack(1);
	}

	@Override
	public IVector getCrackOpeningParameter() {
		return Xi;
	}

	@Override
	protected double computeDeterminantOfCrackJacobian(double... s) {
		return crackLength / 2;
	}

	/**
	 * compute J operator
	 * 
	 * @param s
	 *            local coordinate on the crack
	 * @return J operator
	 */
	protected IMatrix computeCrackInterpolationMatrix(double s) {
		int dim = getMeshPart().getDimension();
		IMatrix J = new Array2DMatrix(dim, dim);
		J.set(0, 0, normal[0]);
		J.set(1, 0, normal[1]);
		J.set(0, 1, tangent[0]);
		J.set(1, 1, tangent[1]);
		return J;
	}

	@Override
	protected void computeEnhancedStrainOperator(MeshPart mp, IMatrix g, double... xi) {
		int dim = mp.getDimension();
		int nn = mp.countModes();
		int ntrs = dim * (dim + 1) / 2;
		Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * nn);
		computeB(mp, B, xi);

		IMatrix Ba = new Array2DMatrix(ntrs, dim);
		for (int i = 0; i < posDomain.size(); i++) {
			int shift = dim * (posDomain.get(i));
			Ba.add(0, 0, -B.get(0, shift));
			Ba.add(1, 1, -B.get(1, 1 + shift));
			Ba.add(2, 0, -B.get(2, shift));
			Ba.add(2, 1, -B.get(2, 1 + shift));
		}
		IMatrix J = computeCrackInterpolationMatrix(0);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, Ba, BLAM.NO_TRANSPOSE, J, 0.0, g);
	}

	@Override
	protected void computeProjectionOperator(MeshPart mp, IMatrix g, double... xi) {
		double coeff = -crackLength / elementArea;
		g.set(0, 0, coeff * (normal[0] * normal[0]));
		g.set(1, 0, coeff * (normal[1] * normal[1]));
		g.set(2, 0, coeff * (2 * normal[0] * normal[1]));
		g.set(0, 1, coeff * (normal[0] * tangent[0]));
		g.set(1, 1, coeff * (normal[1] * tangent[1]));
		g.set(2, 1, coeff * (normal[0] * tangent[1] + normal[1] * tangent[0]));
	}

	@Override
	protected IVector computeCrackTraction(double... s) {
		IVector Tgamma = getCohesiveLaw().computeJumpTraction(this, Xi);
		IMatrix J = computeCrackInterpolationMatrix(0);
		IVector T = new ArrayVector(Xi.getSize());
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, J, Tgamma, 0.0, T);
		return T;
	}

	@Override
	protected IMatrix computeCrackStiffness(double... s) {
		IMatrix Cgamma = getCohesiveLaw().computeJumpStiffness(this, Xi);
		IMatrix J = computeCrackInterpolationMatrix(0);
		IMatrix C = new Array2DMatrix(Xi.getSize(), Xi.getSize());

		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, J, BLAM.NO_TRANSPOSE, Cgamma, 0.0, C);

		// Debugger.watch("J = ", J);
		// Debugger.watch("Cgamma = ", Cgamma);

		// MyLinAlg.btCB(J, Cgamma, 1.0, C);
		// MyLinAlg.bCBt(J, Cgamma, 1.0, C);

		// Debugger.watch("C = ", C);

		return C;
	}

	/**
	 * compute an equivalent stress on crack surface for checking onset criteria
	 * 
	 * @return
	 */
	protected IVector computeEquivalentStressOnCrack() {
		IVector sigma = (IVector) getStress(0).clone();
		int ng = getNumIntegrationPoints();
		for (int i = 1; i < ng; i++) {
			BLAM.add(1.0, getStress(i), sigma);
		}
		BLAM.scale(1.0 / ng, sigma);
		return sigma;
	}

	@Override
	protected boolean checkOnsetCriteria() {
		IVector sigma = computeEquivalentStressOnCrack();

		double oxx = sigma.get(0);
		double oyy = sigma.get(1);
		double oxy = sigma.get(2);
		double t1 = (oxx + oyy) / 2;
		double t2 = Math.sqrt(Math.pow((oxx - oyy) / 2, 2) + Math.pow(oxy, 2));
		double o1 = t1 + t2;

		if (o1 < getCohesiveLaw().getTensileStrength()) {
			return false;
		} else {
			/*
			 * initialize crack related parameters
			 */
			posDomain = new ArrayList<Integer>();
			negDomain = new ArrayList<Integer>();
			crackBoundary = new ArrayList<double[]>();
			Xi = new ArrayVector(2);

			MeshPart mp = getMeshPart();

			/*
			 * compute crack center
			 */
			double[] crackCenter = mp.computeCenter();

			/*
			 * realize crack normal
			 */
			double n1, n2;
			if (Math.abs(oxy) < MathUtilities.tol14) {
				if (oxx > oyy) {
					n1 = 1;
					n2 = 0;
				} else {
					n1 = 0;
					n2 = 1;
				}
			} else {
				double g = Math.sqrt(Math.pow(oxx - o1, 2) + Math.pow(oxy, 2));
				n1 = -oxy / g;
				n2 = (oxx - o1) / g;
			}
			normal = new double[] { n1, n2 };
			tangent = new double[] { -normal[1], normal[0] };

			/*
			 * realize postive domain
			 */
			Node[] nodes = mp.getNodes();
			int nn = mp.countModes();
			// seggregate nodes on one side of the crack
			double a = normal[0];
			double b = normal[1];
			double c = -(a * crackCenter[0] + b * crackCenter[1]);
			for (int i = 0; i < nn; i++) {
				double[] p = nodes[i].getPosition();
				if (a * p[0] + b * p[1] + c > 0) {
					posDomain.add(i);
				} else {
					negDomain.add(i);
				}
			}
			// reverse if PosDomain is not in the direction of normal
			double[] p = nodes[posDomain.get(0)].getPosition();
			if (normal[0] * (p[0] - crackCenter[0]) + normal[1] * (p[1] - crackCenter[1]) < 0) {
				@SuppressWarnings("unchecked")
				ArrayList<Integer> tmp = (ArrayList<Integer>) posDomain.clone();
				posDomain = negDomain;
				negDomain = tmp;
			}

			/*
			 * compute crack length and two intersections
			 */
			// realize two edges intersect with crack
			int[] edge1, edge2;
			if (nn == 3) {
				if (posDomain.size() == 1) {
					edge1 = new int[] { posDomain.get(0), negDomain.get(0) };
					edge2 = new int[] { posDomain.get(0), negDomain.get(1) };
				} else if (posDomain.size() == 2) {
					edge1 = new int[] { negDomain.get(0), posDomain.get(0) };
					edge2 = new int[] { negDomain.get(0), posDomain.get(1) };
				} else {
					throw new Error("Something wrong with the domain");
				}
			} else if (nn == 4) {
				if (posDomain.size() == 1) {
					int cnt = 0;
					int[] adjNodes = new int[2];
					for (int i = 0; i < 3; i++) {
						if (Math.abs(posDomain.get(0) - negDomain.get(i)) != 2) {
							adjNodes[cnt++] = negDomain.get(i);
						}
					}
					edge1 = new int[] { posDomain.get(0), adjNodes[0] };
					edge2 = new int[] { posDomain.get(0), adjNodes[1] };
				} else if (posDomain.size() == 2) {
					if (Math.abs(posDomain.get(0) - negDomain.get(0)) != 2) {
						edge1 = new int[] { posDomain.get(0), negDomain.get(0) };
						edge2 = new int[] { posDomain.get(1), negDomain.get(1) };
					} else {
						edge1 = new int[] { posDomain.get(0), negDomain.get(1) };
						edge2 = new int[] { posDomain.get(1), negDomain.get(0) };
					}
				} else if (posDomain.size() == 3) {
					int cnt = 0;
					int[] adjNodes = new int[2];
					for (int i = 0; i < 3; i++) {
						if (Math.abs(negDomain.get(0) - posDomain.get(i)) != 2) {
							adjNodes[cnt++] = posDomain.get(i);
						}
					}
					edge1 = new int[] { negDomain.get(0), adjNodes[0] };
					edge2 = new int[] { negDomain.get(0), adjNodes[1] };
				} else {
					throw new Error("Something wrong with the domain");
				}
			} else {
				throw new Error("invalid number of nodes");
			}
			// compute intersection with edge 1
			double[] N1 = nodes[edge1[0]].getPosition();
			double[] N2 = nodes[edge1[1]].getPosition();
			double xi1 = N1[0], yi1 = N1[1];
			double xi2 = N2[0], yi2 = N2[1];
			double a1 = yi2 - yi1;
			double b1 = xi1 - xi2;
			double c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
			double[] X1 = MathUtilities.computeIntersection(a, b, c, a1, b1, c1);
			crackBoundary.add(X1);
			// compute intersection with edge 2
			N1 = nodes[edge2[0]].getPosition();
			N2 = nodes[edge2[1]].getPosition();
			xi1 = N1[0];
			yi1 = N1[1];
			xi2 = N2[0];
			yi2 = N2[1];
			a1 = yi2 - yi1;
			b1 = xi1 - xi2;
			c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
			double[] X2 = MathUtilities.computeIntersection(a, b, c, a1, b1, c1);
			crackBoundary.add(X2);
			// compute crack length
			crackLength = LinAlg.diffNorm2(2, X1, X2);

			/*
			 * compute area
			 */
			elementArea = getFace().computeArea();

			return true;
		}
	}
}
